% This program conducts the predictive regression analysis.
%
% 10-10-07 by Shu Yan.

clear;close all;pack;

load ml_data_5                                    % Results for implied vol errors are in ml_04_YZ_022605

R=ret;
n=length(ret);                                    % length of the sample

dt=1/52;                                          % Delta t
tau=1/12;                                         % Investment horizon
n=length(R);

Parm=zeros(12,3);                                 % Initialize the parameter estimates

% The estimates of "ml_102_4.m" - SV-SJ model.
load ml_102_5 parm fval hessian
m_Q=parm(1);s_Q=parm(2);
load ml_102_4_YZ Y0 Z0 parm0
gam=1.91744;
Parm(:,3)=[parm0(1:2) parm0(5) parm0(3:4) parm0(6) m_Q s_Q parm0(7:9) gam]';
Y_svsj=Y0;Z_svsj=Z0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Construct the equity premium %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We only need to find B and C since the parameters are under the objective probability measure.
m_Y=Parm(1,3);k_Y=Parm(2,3);s_Y=Parm(3,3);m_Z=Parm(4,3);k_Z=Parm(5,3);s_Z=Parm(6,3);m_Q=Parm(7,3);s_Q=Parm(8,3);r_12=Parm(9,3);r_13=Parm(10,3);r_23=Parm(11,3);
Pi=[m_Y;m_Z];Lambda=[k_Y, 0;0, k_Z];Gamma=[s_Y^2,r_23*s_Y*s_Z;r_23*s_Y*s_Z,s_Z^2];
Theta=[-0.5*gam*(gam-1),0;
        0,exp(-gam*log(1+m_Q)+0.5*gam*(gam-1)*s_Q^2)*((gam*(1+m_Q)-(gam-1)*exp(gam*s_Q^2))-1)];
[T,Y]=ode23(@ml_ode_obj,[0 tau],[0 0 0 0 0],[],Pi,Theta,Lambda,Gamma);
m=length(T);
B=Y(m,1:2)';C=[Y(m,3) Y(m,4);Y(m,4) Y(m,5)];

prem1=gam*Y0.^2-[r_12*s_Y,r_13*s_Z]*B*Y0-2*[r_12*s_Y,r_13*s_Z]*C(1,:)'*Y0.^2;
prem2=-2*[r_12*s_Y,r_13*s_Z]*C(2,:)'*Y0.*Z0;
prem3=-Z0.^2.*(exp(-gam*log(1+m_Q)+0.5*gam*(gam-1)*s_Q^2)*(1+m_Q-exp(gam*s_Q^2))-m_Q);
prem=prem1+prem2+prem3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Regression Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Regout=zeros(7,4);
% We first consider regressing returns on lagged state variables and equity premium.
Beta=zeros(3,4);Tstat=Beta;Rsqr=zeros(4,1);
y=(exp(ret)-1)*52;
y4=filter(ones(4,1)/4,1,y);                      % Moving average of returns
y8=filter(ones(8,1)/8,1,y);
y13=filter(ones(13,1)/13,1,y);
for I1=1:4
    if I1==1
        yy=y(2:n);x=[ones(n-1,1),Y_svsj(1:n-1)];betaid=[1];
    elseif I1==2
        yy=y(2:n);x=[ones(n-1,1),Z_svsj(1:n-1)];betaid=[2];
    elseif I1==3
        yy=y(2:n);x=[ones(n-1,1),Y_svsj(1:n-1),Z_svsj(1:n-1)];betaid=[1,2];
    elseif I1==4
        yy=y(2:n);x=[ones(n-1,1),prem(1:n-1)];betaid=[3];
    end
    nlag=1;rst=nwest(yy,x,nlag);
    Beta(betaid,I1)=rst.beta(2:(size(x,2)));
    Tstat(betaid,I1)=rst.tstat(2:(size(x,2)));
    Rsqr(I1)=rst.rsqr;
end

Regout(1:2:5,1:4)=Beta;
Regout(2:2:6,1:4)=Tstat;
Regout(7,1:4)=Rsqr';

info.fmt='%3.3f';
info.rnames=strvcat(' ','$Y$ &',' &','$Z$ &',' &','$\phi$ &',' &','$R^2$ &');
disp(['Regression results for lag = ' num2str(nlag)]);
lprint(Regout,info);

% We next consider other macro-variables that are related to stock returns.
dat=dlmread('data2.txt',','); %[divy;pe;tb3m;tb1y;tb10y;aaa;bbb;swap10ymid;swap10yoff;libor3m;cpf3m;cpn3m]
dat=dat(:,1:366)';
dat(find(dat==-999))=NaN;
div=dat(:,1)/100;                              % dividend yield
pe=log(dat(:,2));                              % log pe ratio
tb3m=dat(:,3)/100;                             % 3-month tbill rate
tmsprd=(dat(:,5)-dat(:,4))/100;                % term spread = 10-year yield - 1-year yield
cdsprd=(dat(:,7)-dat(:,6))/100;                % credit spread = bbb yield - aaa yield --- only available from end of 96
swsprd=(dat(:,8)-dat(:,5))/100;                % swap spread = 10-year swap rate - 10-year yield
tdsprd=(dat(:,10)-dat(:,3))/100;               % ted spread = 3-month LIBOR rate - 3-month tbill rate
cfsprd=(dat(:,11)-dat(:,3))/100;               % comercial paper financial spread = 3-month cpf rate - 3-month tbill rate
datmac=[tb3m,cfsprd];
datstats=[];
for i=1:2
    y=datmac(:,i);if i==2,y=y(54:366);end
    auto_y=autocorr(y);
    datstats=[datstats,[mean(y);std(y);skewness(y);kurtosis(y);max(y);min(y);auto_y(2)]];
end
info.fmt='%3.3f';
info.cnames=strvcat('& Mean','Std.','Skew.','Kurt.','Max.','Min.','Autocorr.');
info.rnames=strvcat(' ','Tbill &','CPS &');
disp(['Stats of macro-variables:']);
lprint(datstats',info);

% Regress changes of phi, and V and Lambda on changes of macro-variables.  
Beta=zeros(2,5);Tstat=Beta;Rsqr=zeros(2,4);
for I0=1:4
    if I0==1
        x=diff(Y_svsj.^2)/52;
    elseif I0==2
        x=diff(Z_svsj.^2)/52;
    elseif I0==3
        x=diff([Y_svsj.^2,Z_svsj.^2])/52;
    else
        x=diff(prem)/52;
    end
    for I1=1:2
        if I1==2
            nlag=4;y=diff(datmac(54:n,I1));rst=nwest(y,[ones(n-54,1),x(54:n-1,:)],nlag);
        else
            nlag=4;y=diff(datmac(1:n,I1));rst=nwest(y,[ones(n-1,1),x],nlag);
        end
        if I0<=2
            Beta(I1,I0)=rst.beta(2);Tstat(I1,I0)=rst.tstat(2);
        elseif I0==3
            Beta(I1,3:4)=rst.beta(2:3)';Tstat(I1,3:4)=rst.tstat(2:3)';
        elseif I0==4
            Beta(I1,5)=rst.beta(2);Tstat(I1,5)=rst.tstat(2);
        end
        Rsqr(I1,I0)=rst.rsqr;
    end
end
Regout=zeros(14,2);
Regout(1,:)=Beta(:,1)';Regout(2,:)=Tstat(:,1)';Regout(3,:)=Rsqr(:,1)';
Regout(4,:)=Beta(:,2)';Regout(5,:)=Tstat(:,2)';Regout(6,:)=Rsqr(:,2)';
Regout(7:2:9,:)=Beta(:,3:4)';Regout(8:2:10,:)=Tstat(:,3:4)';Regout(11,:)=Rsqr(:,3)';
Regout(12,:)=Beta(:,5)';Regout(13,:)=Tstat(:,5)';Regout(14,:)=Rsqr(:,4)';
info.fmt='%3.3f';
info.rnames=strvcat(' ','$\Delta V$ &',' &','$R^2$ &','$\Delta\lambda$ &',' &','$R^2$ &','$\Delta V$ &',' &','$\Delta\lambda$ &',' &','$R^2$ &','$\Delta\phi$ &',' &','$R^2$ &');
info.cnames=strvcat('& $\Delta$Tbill','$\Delta$CPS');
disp(['Regress changes of phi, V and Lambda on changes of macro-variables:']);
lprint(Regout,info);

